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ABSTRACT 

We^begin  the  construction  and  the  analysis  of  nonoscillatory  shock  capturing 
methods  for  the  approximation  of  hyperbolic  conservation  laws.  These  schemes 
share  many  desirable  properties  with  total  variation  diminishing  schemes,  but 
TVD  schemes  have  at  most  first  order  accuracy,  in  the  sense  of  truncation  error, 

s  l" 

at  extrema  of  the  solution.  In  this  paper  we  construct:^  uniformly  second  order 
approximation,  which  is  nonoscillatory  in  the  sense  that  the  number  of  extrema  of 
the  discrete  solution  is  not  increasing  in  time.  This  is  achieved  via  a  non- 
oscillatory  piecewise  linear  reconstruction  of  the  solution  from  its  cell  averages, 
time  evolution  through  an  approximate  solution  of  the  resulting  initial  value 
problem,  and  averaging  of  this  approximate  solution  over  each  cell. 


m W.*  k _ L 


'»v>  i  _J  X'fc  j  '~\j.  ~ 


1 


<  '  f  '  t 

AMS  (MOS)  Subject  Classifications:  Primary  65M10,  Secondary  65M05 

Key  Words:  Conservation  Laws,  Finite  Difference  Scheme,  Nonoscillatory,  TVD. 

Work  Unit  Number  3  -  Numerical  Analysis  and  Scientific  Computing 


^School  of  Mathematical  Sciences,  Tel-Aviv  University,  Israel  and  Department  of 
Mathematics,  University  of  California,  Los  Angeles,  CA  90024. 

2 

Department  of  Mathematics,  University  of  California,  Los  Angeles,  CA  90024. 

1)  Research  sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-80-C-0041 
while  in  residence  at  the  MRC,  University  of  Wisconsin-Madison,  and  by  NASA 
Consortium  Agreement  INCA2-IR390-403  and  ARO  Grant  IDAAG29-82-0090  while  at  UCLA. 

2)  Research  supported  by  NSF  Grant  #82-00788,  ARO  Grant  #DAAG29-82-0090,  NASA 
Consortium  Agreement  #NCA2-IR390-403,  and  NASA  Langley  Grants  #NAG  1-270,  and 
NAGI-506. 


UNIFORMLY  HIGH-ORDER  ACCURATE  NON-OSCILLATORY  SCHEMES  I, 

1  2 
Ami  Harten  and  Stanley  Osher 


1.  Introduction.  In  this  paper  we  consider  numerical  approximations  to  weak 
solutions  of  the  scalar  initial  value  problem  (IVP) 


(1.1a) 


U,  +  f(u)s  =  U,  +  a(u)  ux  =  0 


(1.1b) 


u(x,0)  =  Uq(x) 


The  initial  data  u0(x)  are  assumed  to  be  piecewise-smooth  functions  that 
are  either  periodic  or  of  compact  support. 

Let  vf  =  vk(xj,  t„),  xj  =  jh,  t„  —  nr,  denote  a  numerical  approximation  in 


conservation  form 


(1.2a) 


v7+1  -  vf  -  *(fj+V 2  ~  fj-V2)  m  (Eh  *  vn)j 


Here  Eh  is  the  numerical  solution  operator,  X  =  r/h,  and  fj+m>  the 


numerical  flux  is  a  function  of  2k  variables 


(1.2b) 


f]+\n  -  /(v"-i+i,...,v;+t) 


which  is  consistent  with  (1.1a)  in  the  sense  that 


(1.2a) 


/(u, «,...,«)  =  f(u )  . 
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We  consider  the  numerical  approximation  vh(x,t)  in  (1.2)  to  be  a1 
piecewise -constant  function 

(1.3)  vh(x,t)  =  vj,  Xj-u2  <x<xJ+v2,  nr  <  t  S  (n  +  1)t  . 

/)-/ 

Accordingly  we  define  its  total  variation  in  x  to  be  u — 
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(1.4) 


TV(y”)  -  7V(v, (•,/„))  =  2  |v/+1  -  v/|  . 

If  the  total  variation  of  the  numerical  solution  is  uniformly  bounded  in  A  for 

(1.5)  7V(v, (•,/))  <;  C  *  'iV(uo) 

then  any  refinement  sequence  A  -*  0,  t  =  0(A)  has  a  subsequence  Ay  -  0  so 
that 

(1.6)  vhj - >  u 

where  u  is  a  weak  solution  of  (1.1). 

If  all  limit  solutions  (1.6)  of  the  numerical  solution  (1.2)  satisfy  an  entropy 
condition  that  implies  uniqueness  of  the  IVP  (1.1),  then  the  numerical  scheme  is 
convergent  (see  e.g  [3],  [12]). 

Recently  we  have  introduced  the  notion  of  total  variation  diminishing  (TVD) 
schemes  (see  [3]),  where  the  approximate  solution  operator  is  required  to  dimin¬ 
ish  the  total  variation  (1.4)  of  the  numerical  solution  at  each  time-step 

(1.7)  7V(v"+1)  s  TV(v")  ; 

these  schemes  trivially  satisfy  (1.5)  with  C  -  1.  Some  early  work  along  these 
lines  was  done  by  van  Leer  in  [15]. 

TVD  schemes  are  non-osdllatory  in  the  sense  that  the  number  of  local 
extrema  in  the  numerical  solution  is  diminishing  in  time  (as  is  customary  we  use 
"diminishing"  loosely  as  short  for  "non-increasing",  throughout  this  paper). 
Moreover,  the  value  of  an  isolated  local  maximum  may  only  decrease  in  time, 
while  that  of  a  local  minimum  may  only  increase. 
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We  were  able  to  construct  TVD  schemes  that  in  the  sense  of  local  truncation 
error  are  high-order  accurate  everywhere  except  at  local  extrema  where  they 
necessarily  degenerate  into  first-order  accuracy  (see  [4],  [13],  [10],  [11],  [14]). 
The  perpetual  Hamping  of  local  extrema  determines  the  cumulative  global  error  of 
the  "high-order  TVD  schemes"  to  be  0(h)  in  the  L«  norm,  0(h3r2)  in  the  L2 
norm  and  0(h2)  in  the  Lj  norm  (see  [17]). 

In  this  paper  we  introduce  a  larger  class  of  non-osdllatory  schemes,  in  which 
the  solution  operator  is  only  required  to  diminish  the  number  of  local  extrema  in 
the  numerical  solution.  Unlike  TVD  schemes,  which  are  a  subset  of  this  class, 
non-osdllatory  schemes  are  not  required  to  damp  the  values  of  each  local 
extremum  at  every  single  time-step,  but  are  allowed  to  occasionally  accentuate  a 
local  extremum. 

In  a  sequence  of  papers,  of  which  the  present  paper  is  the  first,  we  show  how 
to  construct  non-osdllatory  schemes  that  are  uniformly  high-order  accurate  (in 
the  sense  of  global  error  for  smooth  solutions  of  (1.1)).  In  this  first  paper  we 
describe  a  second-order  accurate  scheme  of  this  type. 

The  fact  that  the  number  of  local  extrema  in  the  numerical  solution  may  only 
diminish  in  time  is  suffident  by  itself  to  guarantee  that  the  application  of  the 
scheme  to  monotone  data  results  in  a  monotone  function.  Thus  non-osdllatory 
schemes,  like  TVD  schemes,  are  monotonidty  preserving.  In  particular,  when 
applied  to  a  step-function,  they  do  not  generate  spurious  oscillations. 

We  note  that  since  the  number  of  local  extrema  in  the  solution  of  non- 
osdllatory  schemes  is  bounded  by  that  of  the  initial  data,  uniform  boundedness  of 
its  total  variation  (1.5)  follows  immediately  if  the  maximum  norm  of  the  solution 
is  shown  to  be  uniformly  bounded. 


2.  Design  Principle  and  Overview 

In  this  section  we  describe  how  to  construct  a  non-osdllatory  scheme  that  is 
uniformly  second-order  accurate. 

Integrating  the  partial  differential  equation  (1.1a)  over  the  computational  cell 
*j+ 1/2)  x  (*n,  *,,+1)  we  get 

(2.1a)  -  S?  -  M/y+wM  -  /,.„(«)] 

where 

(21b)  fj+vM  =  -■  f(u(xj+v2,  0)  dt , 

and 

(2*lc)  J  • 

We  observe  that  although  (2.1a)  is  a  relation  between  the  cell-averages  uij 
and  uy+1,  the  evaluation  of  the  fluxes  fj+y 2(h)  in  (2.1b)  requires  knowledge 
of  the  solution  itself  and  not  its  cell-averages. 

As  in  Godunov's  scheme  and  its  second-order  extension  by  van  Leer  [16]  and 
Colella  and  Woodward  [2],  we  derive  our  scheme  as  a  direct  approximation  to 
(2.1).  We  denote  by  vf  the  numerical  approximation  to  the  cell-averages  uf  of 
the  exact  solution  in  (2.1c),  and  set  vf  to  be  the  cell-averages  of  the  initial  data. 
Given  v"  =  {vf}  we  compute  v*+1  as  follows: 

First  we  reconstruct  u(x,t„)  out  of  its  approximate  cell-averages  {vf}  to  the 
appropriate  accuracy  and  denote  the  result  by  L(x;  v”).  Next  we  solve  the  IVP 

(2  2)  v,  +  f(v)x  =  0,  v(*,0)  =  L(x;  v") 
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and  denote  its  solution  by  v(x,r).  Finally  we  obtain  vjl+1  by  taking  cell- 
averages  of  v(x,t) 

(2.3)  **'-}£? 

The  averaging  operator  in  (2.3)  is  non-osdllatory,  therefore  the  number  of 
local  extrema  in  vn+1  (interpreted  as  a  mesh-function  or  the  piecewise-constant 
function  (1.3))  does  not  exceed  that  of  v(x,t).  Assuming  v(x,f)  to  be  the  exact 
solution  of  (2.2)  implies  (since  the  exact  solution  operator  is  TVD)  that  the 
number  of  local  extrema  in  v(x,t)  is  less  than  or  equal  to  that  of 
*(x,0)  =  L(x;  v").  Therefore  if  the  number  of  local  extrema  in  L(x;  vn)  does 
not  exceed  that  of  v*,  then  the  resulting  scheme  is  non-osdllatory. 

We  condude  that  the  design  of  non-osdllatory  high  order  accurate  schemes 
essentially  boils  down  to  a  problem  oa  the  level  of  approximation  of  functions: 
Given  cell-averages  Uj  of  a  piecewise-smodth  function  u(x),  reconstruct  u(x) 
to  a  desired  accuracy.  Prior  to  studying  this  problem  we  tackle  another  related 
question  in  approximation  of  functions,  that  of  constructing  a  non-osdllatory 
high-order  accurate  interpolation  of  piecewise-smooth  functions. 

In  section  3  we  construct  a  non-osdllatory  piecewise-  parabolic  function 
Q(x;  u )  that  interpolates  a  piecewise-smooth  function  u(x)  at  the  mesh  points 

(2.4a)  Q(xy,  u )  =  u(xj) 

and  satisfies,  wherever  u(x)  is  smooth, 

(2.4b)  Q(x;  u)  =  u(x)  +  0(h*) 

(2-4c)  ~dx  ^ X  ±  °’  =  ~dx  +  ' 
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In  section  4  we  make  use  of  this  non-osdUatory  piecewise- parabolic  interpo¬ 
late  to  design  a  non-osdUatory  reconstruction  of  a  piecewise-smooth  function 
from  its  ceU-averages.  As  in  [16],  [2],  [5],  and  [9]  we  take  L(x:  u)  to  be  the  fol¬ 
lowing  piecewise-linear  function 

(2.5a)  L(x;  u)  =  uj  +  Sj(x  ~  xj)/h  for  \x  -  xj\<  h/2  . 

Unlike  the  above  references  that  present  "second-order  accurate"  TVD 
schemes,  we  compute  the  slopes  SJh  from  Q(x;  it)  by 

(2.5b)  SjJh  =  m|-£-  Q(xj  -  0;u),  Q(xj  +  0;  u)j  - 

(fere  m(x,y)  is  the  min  mod  function 

(s  •  min(|x|,iyD  i*  *gn(*)  =  **“00  *  s 

(2.6)  m(x,y)  ~  iq  otherwise 

We  show  in  section  4  that  L(x;u)  is  a  proper  reconstruction  of  u(x)  in  the 
sense  that  whenever  u(jc)  is  smooth 

(2.7a)  L(x;  u)  =  «(x)  +  O^h2) 

and 

(2.7b)  L(x;  5)  =  u(x)  +  0{h?)  . 

Here  u(x)  =  h'1  u(x  +  y)  dy  and  L(x;  u)  - 
=  h'1  L(x  +  y;  u)  dy  ;  like  Q(x;  u),  the  latter  is  also  a  non-osdllatory 
piecewise-parabotic  interpolant  of  u(x), 

(2.7c)  L(xj;  u)  =  u(xj)  . 
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We  remark  that  the  "second-order  accurate"  TVD  schemes  described  in  the 
above  mentioned  references  use  a  slope  Sj/h  in  (2.5a)  that  approximates 
(d/dx)  u(xj )  to  0(h),  and  their  loss  of  second-order  accuracy  at  local  extrema 
points  is  due  to  lack  of  smoothness  of  the  coefficient  in  the  0(h)  term  at  these 
points  *.  This  problem  is  circumvented  in  the  present  scheme  by  taking  Sj/h  to 
be  (2.5b)  which  is  an  0(h?)  approximation  to  (d/dx)  u(xj).  Unfortunately  there 
is  a  price  to  pay  for  this  extra  accuracy,  namely  the  loss  of  the  TVD  property.  As 
in  TVD  schemes 

(2.8)  TV(vn")  j£  TV(L(",  v"))  , 

however  here 

TV(L(  \  v"))  te  TV(vn) 

and  indeed  the  scheme  may  occasionally  increase  the  variation  of  the  numerical 
solution.  Although  we  prove  that  the  scheme  is  non-osdllatory  we  have  not  been 
able  as  yet  to  complete  a  proof  of  uniform  boundedness  of  the  total  variation  of 
the  numerical  solution;  this  is  due  to  lack  of  techniques  to  verify  uniform  bound¬ 
edness  of  the  maximum  norm  of  the  numerical  solution. 

In  section  5  we  study  the  proposed  scheme  in  the  constant  coefficient  case. 
We  verify  that  it  is  uniformly  second-order  accurate,  examine  its  behavior  at  local 
extrema  points  and  get  estimates  for  the  possible  increase  in  total  variation  per 
time-step. 

In  this  paper  where  we  consider  numerical  schemes  of  the  form  (1.2)  that  are 

1.  We  repeat  that  the  results  of  [8]  and  [11]  imply  that  TVD  schemes,  no  matter  how 
they  are  constructed,  must  have  this  loss  of  accuracy  at  local  extrema 
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derived  from  approximating  the  relation  (2.1),  it  is  only  natural  to  consider  trun¬ 
cation  error  in  the  sense  of  cell-averages,  i.e.  we  say  that  the  scheme  (1.2)  is 
second-order  accurate  if 

(2.9)  iT*1  =  Ek  ■  uT  +  0(h 3) 

where  u"  is  the  cell-average  (2.1c)  of  the  exact  solution.  Since 

(2.10)  u(x)  =  u(x)  +  Oik2) 

whenever  u(x)  is  smooth,  (2.9)  holds  also  for  pointwise  values  of  the  solution. 
However,  in  the  context  of  3rd  and  higher  order  accurate  schemes,  this  difference 
in  definitions  of  truncation  error  will  be  not  only  conceptual  but  of  practical 
importance  as  well. 

Up  to  this  point  we  have  assumed  that  v(x,t)  in  (2.3)  is  the  exact  solution 
to  (2.2).  The  resulting  scheme 

(2.11a)  M /y+wM  -  fj-v&y) , 

where  fj^m(v)  is  (2-lb)  applied  to  v(x,/), 

(2.11b)  /j.i«(v)  -  i  £  /(v(*.i))  * . 

is  certainly  second-order  accurate  in  the  sense  of  (2.9).  Starting  with  the  exact 
cell- averages  vj  =  Sj  in  (2.11)  we  get  from  (2.7a)  that 

(2.12a)  v(jr,r)  =  u(x,t  +  t„ )  +  0(k*)  for  Osjst 

and  consequently 

(2- 12b)  f)+\dy)  =  fj+v 2(«)  +  0(h*)  , 

which  implies  (2.9)  due  to  the  sufficient  smoothness  of  the  coefficient  in  the 
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Oih2)  term  in  (2.12b). 

In  section  6  we  replace  the  exact  solution  v(x,t)  in  (2.3)  by  an  approximate 
one,  which  we  denote  by  vn(x,t)  .  This  approximate  solution  is  conservative, 
TVD,  and  second-order  accurate  in  the  sense  of  (2.12a).  Thus  replacing  v(x,t) 
in  (2.3)  by  this  approximate  solution  results  in  a  conservative  scheme  that  is  non- 
osdllatory  and  uniformly  second  order  accurate. 

We  remark  that  an  alternative  approach  to  the  above  is  to  approximate 
fj+V2(y)  in  (2.11b)  by  using  a  midpoint  rule  (or  trapezoidal  rule)  for  the 
integral  and  by  replacing  v(x,t)  with  a  non-oscillating  second-order  accurate 
approximate  one  v„(jc,r)  (see  [16]  and  [2]).  The  resulting  scheme 

(2.13a)  v?+1  =  v/  -  -  fj-m) 

(2-  13b)  fj+tn  =  fM*j+ 1/2.  t/2)) 

is  certainly  second-order  accurate,  and  it  is'non-oscillatory  in  the  constant  coeffi¬ 
cient  case.  Since  we  have  not  used  the  cell-averaging  (2.3)  to  derive  this  scheme, 
we  cannot  ascertain  in  general  that  the  resulting  scheme  is  non-osdllatory. 
Nevertheless,  our  numerical  experiments  as  well  as  many  other  experiments  in 
the  context  of  TVD  schemes  (see  e.g.  [1],  [2])  demonstrate  that  the  numerical 
results  are  non-osdllatory  in  many  (if  not  all)  applications. 

In  section  7  we  present  some  numerical  experiments  that  compare  the  present 
scheme  with  a  typical  "second-order  accurate"  TVD  scheme. 

3.  Nonoscillatory  interpolation. 

The  oscillatory  nature  of  second  order  accurate  Lax-Wendroff  type  schemes 
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results  from  a  Gibbs  phenomenon  associated  with  high-order  interpolation  across 
discontinuities.  In  this  section,  as  a  preparatory  step  towards  designing  a  non  os¬ 
cillatory  approximation  to  (1.1.),  we  construct  a  non-oscillatory  piece  wise- 
parabolic  interpolant  Q(x;  u)  to  a  piece  wise-smooth  function  u(x)  such  that 

(3.1a)  Q(xt;  u)  =  u(x,) 

(3.1b)  Q(x\  u)  m  qi+v2(x;  u),  S  i  S  xl+1  , 

where  qi+in  is  a  quadratic  polynomial,  and 

(3.1c)  Q(x;  u )  -  u(x)  =  0(h3),  Q(x  ±  0;  u)  -  ■—  u(x)  =  0(h*) 
wherever  u(x)  is  smooth. 

Q(x;  u )  is  non-oscillatory  in  the  sense  that  the  number  of  its  local  extrema 
does  not  exceed  that  of  u(x). 

Since 

«)  =  “i>  ?i+i/2(*i+ iJ  «)  = 

it  can  be  written  in  the  form 

(3.2a )qi+v2(xvt)  =  «,  +  di+in  u(x  ~  xd/fl  +  yA+i/2  “  ( x  -  x,)(x  -  xl+1)/A2 
where 

(3.2b)  di+y2  u  =  u)+1  -  u, 

and  Di+y2u  is  yet  to  be  determined. 

A +i/2  u  =  h2  q-+ i/2(x;  u),  xf  s  x  s  xf+1  . 
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Q  has  its  local  extremum  in  the  interior  of  some  (x„  x,+ j).  We  recall  that  such 
an  extremum  is  characterized  by  ki+l/2  v"j  <  1^2  |D,+  i/2  v#*I  and  that  Sf  and 
Sf+1  in  this  case  are  given  by  (4.14);  therefore  Sf+i  -  Sf  =  Df+1/2  v".  From 
(5.3a)  and  (4.6a)  we  see  that  in  general 

(5. 10a)  vf^1  -  Q(x.+i  -  at;  v")  =  y  n(l  -  p.)(A+y 2  v"  -  Sf+1  +  Sf)  - 
Hence 

(5.10b)  K+V2  A<\  \Di+in  v"|  =>  vf+"11  =  Q(*i+i  ~  v")  . 

Relation  (5.10b)  confirms  the  second  order  accuracy  of  the  scheme  at  local 
extrema.  Although  it  does  not  necessitate  accentuation  of  the  extremal  values,  as 
vf#  in  (5.10b)  may  still  be  in  [vf,  vf+1],  it  does  allow  vf+i1  to  deviate  from 
this  interval  by  as  much  as 

(5.10c)  i  |Di+w  V|(to+  M  *"m+vi  A  -  Vtf 

Thus  (5.10b)  is  the  essential  difference  between  the  present  scheme  and  the 
“second  order"  TVD  schemes. 

A  similar  analysis,  which  we  do  not  present  here,  shows  that  if  vf  is  a 
mesh-extremum  then  vf*1,  j  =  «',  i  +  1,  relates  to  Q(xj  -  ar;  v")  in  the  fol¬ 
lowing  way: 

(5.11a)  vf*1  ^  Q(xj  -  at;  v"),  j  =  i,  i  +  1,  if  vf  is  a  maximum 
(5.11b)  v"+1  s  Q(xj  —  at;  v"),  j  =  i,  i  +  1,  if  vf  is  a  minimum  . 

From  (5.9)-(5.11)  we  deduce  the  following  relation  between  the  total  varia- 


-24- 


(5.6b)  \ij-v a  v"|  as  |  1 Dj-v2  v*|  =>  by-1/2 I  -  1*7  "  V|  *  2 

Rewriting  (5.3)  in  this  case  as 

(5.7a)  v;+1  -  v/  -  n  dy-i/2  v*  -  1/2  n(l  -  »)yj-v2dj-m  vH 
m  (1  -  <*y-l/2)  v"  +  v"-l 

where 

(5.7b)  =  M-  +  y  P-C1  ~  P->Yy-l/2  • 

We  see  that  the  CFL  condition  0  <  u.  ^  1  and  (5.6b)  imply  that 
(5.7c)  0  s  <r;_i/2  ^  1  ; 

thus  we  conclude 

(5.8)  \ij-ut  v“|  z  i  pj-m  ."I  =»  »7+1  €  [ij-i,  vfl . 


Relation  (5.8)  shows  that  if  v"  is  monotone  for  JLi£  j  &  JR,  i.e. 

VJL  s  vytM  2:  •  •  •  2S  V7jr,  or  Vy4  a  Vy^,  a  •  •  •  ^  vjM,  then  vfl+1  is  mono¬ 
tone  for  4  +  1  s  j  s  JR,  and  in  the  same  sense.  Relation  (5.8)  also  shows 
that  mesh-extrema  of  v",  i.e.  those  for  which  Q  has  its  local  extremum  at  a 
mesh  point,  arc  being  damped  at  the  n-th  time-step.  Namely, 


(5.9a)  | <tiiln  v"|  2  y  pi±u2  *1,  vj-t  s  v; 


raax(v"  *  . 


v;*1) 


(5.9b)  yjit/1y- 1  a  j  1D/*m  n 


We  turn  now  to  consider  interior  local  extrema  of  v",  i.e.  those  for  which 


(5.5a)  7V(vn+1)  s:  7V(L(*;  v"))  . 

Using  (4.15)  and  (5.5a)  we  get  the  following  upper  bound  on  the  possible  growth 
of  the  total  variation  of  the  numerical  solution  per  time-step 

(5.5b)  7V(v"+1)  -  TV(v")  fi  2  (i  »"l  -  K.+1/1  «"l]  • 

Here  Mn  is  the  set  of  indices  of  intervals  (xm,  xm+1)  in  the  interior  of 
which  L(x\  v")  (and  also  Q(x;  v1*))  has  a  local  extremum.  The  number  of  these 
intervals  is  finite  and  remains  uniformly  bounded  in  time  by  the  number  of  l<yal 
extrema  in  the  initial  data. 

Clearly  the  upper  bound  (5.5b)  is  overly  pessimistic.  It  estimates  the  possi¬ 
ble  increase  in  variation  in  the  reconstruction  step  due  to  replacing  the  cell- 
averages  vj  by  the  piecewise-linear  function  L(x;  v").  It  does  not  take  into 
account  the  possible  decrease  in  variation  in  the  averaging  step  (2.3),  resulting 
from  doing  just  the  opposite,  i.e.  replacing  the  piecewise-linear  function 
L(x  -  ar;  v")  in  (5.2)  by  its  cell-averages  (5.3a). 

In  the  following  we  shall  examine  the  temporal  behaviour  of  the  local 
extrema  of  the  numerical  solution  and  its  total  variation  by  analysing  the  explicit 
values  of  the  cell-averages  v"+ 5  given  by  (5.3b).  To  simplify  our  presentation 
let  us  assume  that  a  >  0. 

First  we  note  that  (4.8b)  implies 

(5.6a)  \Sj  -  Sj-i |  s  |Sy|  +  \Sj+  ,|  s  2  maxflrf,-^  v"|,  j  | . 

Hence 
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where  L  is  (4.5a).  We  have  shown  in  section  3  that  the  number  of  local 
extrema  in  Lfov*)  does  not  exceed  that  of  v”.  Since  v”*1  in  (5.3a)  is  a  cell- 
average  of  L,  it  follows  that  the  number  of  local  extrema  in  vn+1  does  not 
exceed  that  of  v*.  and  consequently  the  scheme  (5.3a)  is  non-osdllatory. 

Using  (4.5b)  in  (5.3a)  we  get  the  following  expression  for  the  scheme 

(5.3b)  »;+1  =  (£„  •  »")y 

* 

tKtj-vi  v"-v2  tt(i  -  10(57  -  sj-i) 

-  vj-  pdj-nn  **  +  1/2  |i(l  +  |i)(5/+i  -  Sf) 

Eh  denotes  the  operator  form  of  the  finite  difference  scheme; 

CFL-number,  is  assumed  to  satisfy 

(5.3c)  (is.)  as  1 . 

We  turn  now  to  prove  that  (5.3)  is  second  order  accurate  in  the  sense  of 
(2.9),  i.e.  if  u"(x)  denotes  the  mean  value  (4.1)  of  u(x,f„)  then 

(5.4a)  ip1  -  (£»  •  iT),  =  0(/,3) . 

To  show  that  we  observe  that  in  the  constant  coefficient  case  (5.1) 
uj+1  =  u"(xy  -  at),  and  by  (5.3a)  (Eh  •  tT)j  =  L(xj  -  ar;  w”).  Hence  the  LHS 
of  (5.4a)  is  nothing  but 

(5.4b)  xT{xj  -  at)  -  L(xj  -  at;  u")  , 

which  is  0(h3)  as  a  direct  consequence  of  (4.bc). 

Next  we  study  the  time-dependence  of  the  total  variation  and  the  maximum 
norm  of  the  numerical  solution  (5.3).  In  section  2  we  have  pointed  out  that 


if  a  >  0 
if  a  <  0  ’ 

U.  =  Xa,  the 
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(4.14b)  Sj+i  =  Sf+i  -  dj+i/i  u  +  y  Dj+y2  u  . 

The  same  analysis  shows  that  (4.14)  holds  also  for  the  case  that  Q(x;  u)  has  a 
local  minimum  in  (xj,  Xj+l).  (4.9b)  follows  immediately  from  (4.14)  and  (4.7a). 

We  note  that  since  L( x;  u)  is  continuous  at  xj 

(4.15)  7V(i(-;  5))  -  X  ;  i))  =  X  | 

=  X  1/2  “1  +  Pm*  1/2  “  l4n  +  1/2  *nj  • 

Here  M  is  the  set  of  indices  of  intervals  (xmtxm^.{)  in  the  interior  of  which  L 
(and  also  Q)  has  a  local  extremum.  The  number  of  these  intervals  is  finite  and 
is  bounded  by  the  number  of  local  extrema  of  u(x).  Comparing  (4.9)  with  (3.8) 
we  note  that 

(4.16)  7v(M;5))  =  mG(;“))- 

5.  The  constant  coefficient  case. 

In  this  section  we  study  the  constant  coefficient  case 

(5.1)  u,  +  aux  —  0,  a  =  const . 

The  exact  solution  of  the  IV?  (2.2)  is 

(5.2)  v(xft)  =  L(x  -  at;  vn )  . 

Hence  our  scheme  (2.3)  is 

(5.3a)  v/+1  =  J L(x  -  at;  v")  dx  =  L(xj  -  at;  v") 
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To  complete  the  verification  of  (4.7b)  we  still  have  to  show  that 
(4.9b)  kfy+i/Z  **l  <  ~2  ^Pj+y 2  ^  ^  *  **))  —  ~2  1  **l  * 

First  we  observe  that 

(4.10)  Sf  -  Sf  -  (^i*f-y2  “  ~  2"  “)  “  0*1- 1/2  “  +  2"  Di-  1/2  “)  • 

-  DjU  -  y(Dj+i/2  “  +  ^i- 1/2  “) 

Since  (3.2d)  implies 

(4.11)  \Di  51  ^5  y  (lDi-i/2  “1  l°<+ 1/2  “1) 

we  conclude  from  (4.10)  that 

(4.12)  (tf  -  Sf)  sgn(D,  u)*0. 

We  turn  now  to  prove  (4.9b).  First  let  us  consider  the  case  that  Q(x\  u ) 
has  a  local  maximum  in  (xj,Xj+ 1),  i.e.  Dj  u  <  0,  DJ+1  u  <  0,  and 

ty+1/2  **1  <  Y  1 Dj+Vl  “!• 

It  follows  from  (4.12)  that 

(4.13a)  Sf  a  Sf  —  dj+yi  u  -  y  Dj+yi  u  >  0 

(4.13b)  0  >  dj+  !/2  u  +  y  Z)y+i/2  “  =  S/+1  ^  Sf+\  ■ 

The  relations  (4.13)  and  the  definitions  (4.8a),  (4.2b)  imply  that 
(4.14a)  Sj  -  Sf  =  dj+iftit  -  —  Dj+yiu 


-19- 


to  u(x)  in  exactly  the  same  sense  as  Q  is  to  the  interpolated  function  (see  sec¬ 
tion  3). 

Next  let  us  denote 

(4.8a)  Sf  =  h  •  -j-  Q(xj  ±  0;  u)  , 

i.e. 

(4.8b)  Sj~  =  dj-xn  «  +  y  Dj-in  s?  m  dj+vi  “  ~  y  Dj^in  “  . 

and  observe  that  (4.2b)  implies  that 

(4.8c)  y(|Sy|  +  j5/+iD  *  2^  ’  +  lm(S/+l»  ^Al)l] 

*  +  fy+ll)  =  \  [l4+l/2  “  -  \  DJ+V 2  *1  +  to+1/2  «  +  7  S|J 

=  max||dy+!/2  *«1,  Y^P/+1/Z  “lj  • 

We  note  that  if  ty+i/2  51  ^  M2|0/+i/2  then 

sgn|rfj+i/2  u±  y  D)+\n  “j  *  sgn(d/+y2  “)  0 

which  in  turn  implies 

sgn(5j)  *  sgn(<fy+  1/2  5)  a:  0,  sgn(S/+1)  •  sgn(dy+i/2  5)  a  0 

It  follows  then  from  (4.8c)  that  the  RHS  of  (4.7a)  is  ldJ+u2  **!•  Thi*  8ljows  that 

(4.9a)  M/+1/2  «1  *  J  Pj+vi  “1  =»  “»  =  Hf+w  «1  • 
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also  at  local  extrema  points)  and  therefore 


Sj+ 1  -  Sj  ~  h2  u(xJ+v£  +  0(h2)  . 

On  the  other  hand  (3.2)  shows  that 

Dj+V2  “  *  *2  ~2  +  °(fc3)  • 

Therefore 

(4.6b)  Sj+ 1  “  -  £/+ 1/2  “  ~  of*3) 

which  shows  that  RHS  of  (4.6a)  is  0(A)3).  Since  (2.1c)  shows  that 

G(x;  5)  -  u(x)  =  Ofh2) 
we  conclude  from  (4.6a)-(4.6b)  that 
(4.6c)  L(x;  u)  -  u(x)  =  0(h 3)  . 

We  turn  now  to  prove  that  L(x;  u)  is  a  non-osdllatory  approximation  to 
u(x);  this  certainly  implies  that  L(x;  u)  is  a  non-osdllatory  approximation  to 
u(x).  We  shall  do  so  by  showing  that  TV^J>Xj^(L(-;u))r  the  total-variation  of 
L(x\  u)  in  [xj,  x;+1],  which  has  the  value 

(4.7.)  7V[vj.j(£(-;i))  -  A((Sy|  +  (Sy^D  +  | rf/+w«  -  y(Sy  +  Jy+,)l  , 

can  also  be  expressed  as 

(4.7b)  S))  =  mM(|d/+1/2  «1,  y|Dy+  in  JD  - 

Then  it  follows  immediately  form  (4.7b),  (4.3a)  and  (3.8)  that  L  is  monotone  in 
[xj,  xy>1]  if  and  only  if  Q  is;  consequently  L  is  a  non-osdllatory  approximation 


and  therefore  it  follows  from  (3.1c)  that 

(4.4b)  JSJ=-£  “W  +  “  &  “(*/)  +  Of*2)  • 

Consequently  the  RHS  of  (4.2a)  can  be  expanded  as 
(4.4c)  L(x\  u)  =  uj  +  (x  -  xj)  u(xj )  +  OQi2) 

=  u(x)  +  tyh2)  for  jx  -  x}\  <  Y  h 
and  thus  (4.3b)  follows. 

Denote  by  L(x\  u)  the  mean  value  of  L(x ;  u)  in  (x  -  fi/2,  x  +  A/2),  i.e. 
(4.5a)  L(x;  u)  =  j  L{y\ u)  dy  . 

Using  (4.2a)  to  evaluate  the  integral  in  (4.5a)  we  find 
(4.5b)  L(x ;  u)  =  uj  +  dj+v2  u 

'  (*  ~  xjVh  +  (l/2)(S/+i  -  Sj)(x  -  xj)(x  -  xj+1)/h2, 

for  Xj+i 

(4.5c)  L(xj\  u)  =  uj. 

(fence  L(x;  u),  like  Q(x;  u),  is  a  piecewise-parabolic  interpolant  of  u(x). 
Comparing  (4.5b)  with  (3.2)  we  find  that  for  xj  s  x  s  xj+i 

(4.6a)  L(x;  S)  -  Q(x;  u)  =  y(.S/+1  -  Sj  -  DJ+V2  u)(x  -  xj){x  -  xJ+l)/h2  . 
From  (4.4b)  we  see  that  S<  =  h  —  u(x{)  +  0(h3)  (Note  that  this  is  true 


4.  Non-osdllatory  Reconstruction. 

Let  u(x)  be  a  piecewise-smooth  function  and  denote  by  u(x)  its  mean  over 
(x  -  ht 2jt  +  A/2),  i.e. 

(4-i>  -  j  £+“  ^  * 

We  denote  u)  -  u(xj)  and  refer  to  these  values  as  cell-averages  of  u(x).  Given 
{uy},  the  task  in  hand  is  to  reconstruct  u(x)  to  0(h?)  in  a  non-osdllatory  way; 
denote  the  approximately  reconstructed  function  by  L(x;  u).  To  achieve  O^h2) 
accuracy  it  is  sufficient  to  consider  L(x\  u)  to  be  a  piecewise-linear  function.  To 
make  I(x;  u)  a  non-osdllatory  approximation  we  use  the  non-osdllatory  piece- 
wise  parabolic  interpolation  Q(x;  5)  to  compute  its  slopes  as  follows: 

(4.2a)  I(x;  u)  =  u}  +  Sj(x  -  xj)/h  for  \x  -  xy|  <  —  h 

(4.2b)  Sj  =  h  ■  Q(xj  -  0;  u),  Q{x}  +  0;  5)j  . 

Here  m  is  the  min  mod  function  (3.3);  </,+i/2  u  and  £><+  ^  “  are  (3.2b)  and 
(3.2d),  respectively. 

We  note  that  L(x;  5)  may  be  discontinuous  at  {xy+1/2}  and  that 
(4.3a)  L(xj‘,  u)  =  uj  . 

To  see  that  wherever  u(x)  is  smooth 
(4.3b)  L(x,u)  -  u(x)  =  Oik*) 

we  observe  that  in  this  case 
(4.4a)  u(x)  =  u(x)  +  0{h*) 
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(3.8c)  o  s  me)  -  2  Vi+ui -I  =  2v  fc.*» «l  (|^r|  -  i)2 

^  T  2  l°«+i/2  “I  • 

*  miM 

The  sum  in  the  RHS  of  (3.8c)  is  taken  over  the  set  of  indices  M  of  intervals 
(xWf  xm+1)  in  which  fim+m  «|  >  2&,+1/2  u\,  i.e.  where  Q  has  interior- 

Next  we  show  that  if  u(x)  is  a  piecewise-smooth  function  of  bounded  varia¬ 
tion,  then 

(3.9)  lim  7V(Q(*; «))  =  7V(u)  . 

k— 0 

We  observe  that  in  this  case  the  number  of  intervals  in  At  is  finite  and  is  uni¬ 
formly  bounded  by  the  number  of  local  extrema  in  m(x).  Hence  (3.9)  will  follow 
if  we  prove  that  Dm+V2  “  -  0  as  ft  -  0  for  all  m  €  M .  To  accomplish  this  we 
show  that  for  h  sufficiently  small  M  does  not  include  intervals  (xm,  x„,+1)  in 
which  u(x)  is  discontinuous.  To  see  that  let  us  examine  the  case  where  u(x) 
has  a  discontinuity  at  x  €  (xf,  xl+1)  .  Qearly  di+1/2  u  approaches  the  size  of 
the  jump  in  u  while  di-y^u  approaches  zero  as  h-  0.  Consequently 

(3.10a)  ]piu/di+y2  Ml  =  |l  ~  di-yju/di+yi  u|  ■*  1  as  h  —  0 


Hence  for  h  sufficiently  small 

(3.10b)  2\4iJ.y2  “I  >  ^Pi  m|  ^  l^t+l/2  “I 

which  implies  l  g  M. 


mesh  function  {uy}  (for  obvious  reasons  the  case  ut  =  uJ+1  is  counted  as  a 
single-extremum) .  The  above  analysis  also  shows  that  interior-extrema  are  iso¬ 
lated,  i.e.  if  Q  has  an  interior-extremum  in  (xit  xi+1),  then  it  is  the  only  local 
extremum  of  Q  in  fo-i,  x1+2). 

We  turn  now  to  examine  the  case  that  Q  has  a  local  extremum  at  a  mesh 
point  xf;  this  will  be  refered  to  as  a  mesh-extremum.  The  above  observation 
that  interior  extrema  are  isolated  excludes  the  possibility  that  Q  has  an  interior- 
extremum  in  either  (xj-i.Xj)  or  (xit  x1+1)  and  consequently  Q  is  monotone  in 
these  intervals.  This  implies  that  u  •  di^y2  u  <  0  and  therefore  uf  is  a 
local  extremum  of  the  mesh  function  {uy}.  This  concludes  the  proof  that  fi(x;  u) 
is  a  non-osdllatory  interpolant  of  u. 

We  next  express  the  non-osdllatory  nature  of  Q  in  terms  of  total  variation. 
If  Pj+in  “I  ^  2|dy+1/2  u|  then  (2.5)  implies  that  Q  is  monotone  in  (xy,  xy+1]. 
Thus 

(3.8a)  Pj+vi  “I  ^  2|dy+y2  m|  =>  TV[XjtXj^(Q)  =  |rfy+i/2  «| . 

K  IO/+1/2  “I  >  2|dy+1/2  “I  then  Q  has  a  local  extremum  in  (xy,  xy+1)  and 

=  \qj+vi  -  uy|  +  |uy+i  ~  q]±\n\  • 

Using  (2.6)  we  get 

(3.8b)  fDy+j/2  u\  >  2\dj+y2  u|  =>  TV[Xj  X]  ^(Q) 


We  condude  that 


mesh  function  {uj},  the  number  of  which  certainly  does  not  exceed  that  of  u(x). 

Q  may  have  a  local  extremum  in  either  the  interior  of  some  interval 
(*,,*1+1)  or  at  a  mesh  point  x,-.  The  first  case,  which  will  be  fcrcd  to  as 
interior-extremum,  occurs  when  there  is  a  point  x*,  xf  <  x*  <  xj+1,  such  that 

~  Q(x*\ u)  =  0,  but  0  . 

Prom  (3.2a)  if  follows  that  Q  has  an  interior-extremum,  in  (x,,  xi+1)  if  and 
only  if 

(3-5)  \Pi+vi  “I  >  2|rf(+i/2  u| . 

ft+K  =  1/2(0*  ^  value  of  the  interior-extremum  is  then 

(3«)  4f+w  =  ■*  -  J  D1+w«  •  -  y)  ; 

-D/+1/2  a  <  0  it  is  a  local  maximum;  if  Di+y 2  u  >  0  it  is  a  local  minimum . 
Since  Di+V^i  *  m(D(u,  Di+1  u),  (3.5)  holds  if  and  only  if 

(3.7a)  u  •  Dj+fU  >  0 

(3.7b)  \Pj  u|  >  2|d<+j/2  “l»  j  ~  t  +  1  • 

This  implies  that  ft+1/2  has  a  local  extremum  in  (x;,  x/+J)  if  and  only  if  both  q( 
and  qi+\  also  have  a  local  extremum  in  (xj,  xj+i)  and  of  the  same  kind.  Siw* 
a  parabola  has  at  most  one  local  extremum,  it  follows  then  that  qt  does  not  have 
a  local  extremum  in  fa-i,  x,)  and  qi+[  does  not  have  one  in  (xj+1,  xi+2)- 
Consequently  Q  is  monotone  in  both  (x,_!,  x( )  and  (xj+1,  x1+2),  but  in  an 
opposite  sense,  i.e.  di-y^u  •  di+y 2  «  <  0;  the  latter  implies  that  u  has  a  local 
extremum  in  [x(,  xJ+ j]  and  that  either  or  is  a  local  extremum  of  the 


We  consider  as  candidates  for  qt+yi  the  two  quadratic  polynomials  q( 
and  qi+u  interpolating  «(x)  at  (x,_lf  xh  x,+1)  and  (xh  xi+u  x/+2),  respec¬ 
tively,  and  choose  qi+1/2  to  be  the  one  that  is  least  oscillatory  in  [x(,  x(  .J. 
Both  qj,  j  =  i  and  j  =  i  +  1,  can  be  written  as  (3.2a)  with  D^yiu  -  Dju 
where 


(3.2c)  DjU  =  dj+ _  dj—yJU  -  uj+i  -  2uj  +  My_1  . 

Since  the  least  oscillatory  of  qi  and  $}+1  can  be  characterized  as  the  one  that 
deviates  the  least  from  the  line  connecting  (xhu()  with  (x(+1,uj+1)  we  choose 
Di+y iu  in  (3.2a)  to  be 

(3.2d)  A+1/2  “  “  w(2)/iif  Di+iu)  , 

where  m(x,y)  is  the  min  mod  function 


(3.3) 


m(x 


■r)  m  {0 


s  •  min(|x|,lyD  if  sgn(x)  =  sgn(y)  =  s 
otherwise. 


If  u(x)  is  smooth  in  [xy_.lr*y+1]t  then  qj  as  a  quadratic  interpolant  of  u 
satisfies 


(3.4)$y(x)  -  ll(x)  =  0(h\  qj(x)  -  ~  «(x)  =  0(k*),  Xy_j  :£  x  *  Xy  +  1  . 

If  D,u  •  Di+iu  2:  0  then  qt+y 2  is  either  qt  or  qi+y  Otherwise  we  set 
Df+ y 2  «  =  0,  but  then  smoothness  of  u  implies  that  Dju  =  0(/»3)  and  conse¬ 
quently  tfi+i/2  “  qj  ~  0(h2)  for  j  =  i,  i  +  1.  Thus  (3.1c)  follows  from  (3.4). 

We  turn  now  to  prove  that  Q(x;  u)  is  a  nonosdllatory  interpolant  of  u, 
i.e.  that  the  number  of  its  local  extrema  does  not  exceed  that  of  u.  We  do  so  by 
showing  a  one-to-one  correspondence  between  local  extrema  of  Q  to  those  of  the 


tion  of  the  numerical  solution  and  that  of  its  piecewise-  parabolic  interpolant  Q : 
(5.12)  7V({( }(xj  -  ar\  v")})  s  7V(v"+1)  *  TV(Q{  \  v"))  . 

The  LHS  of  (5.12)  is  the  total  variation  of  the  mesh  function 
{Q(xj  -  ot;  v")}.  Relation  (5.12)  suggests  to  consider  an  equivalent  definition 
7V  of  the  total  variation  of  the  numerical  approximation  of  the  form 

fv(yH)  =  TV(G(;v")) 

with  the  hope  that  the  scheme  (5.3)  is  TVD  with  respect  to  this  modified  defini¬ 
tion.  Unfortunately  our  numerical  experiments  have  shown  that  there  are 
instances,  although  rather  rare,  that  TV(yH)  is  increasing  with  it;  the  same  is 
true  for  TV(vn)  =  7V(L(  ;  v")). 

As  we  have  mentioned  in  the  introduction,  because  of  the  nonosdllatory 
nature  of  the  scheme,  uniform  total  variation  boundedness  of  the  numerical  solu¬ 
tion  is  implied  by  uniform  boundedness  of  its  maTimnm  norm.  If  we  follow  a 
particular  local  maximum  of  the  initial  data  we  see  from  (5.9)-(5.10a)  that  it 
actually  decreases  most  of  the  time,  and  whenever  it  does  increase  (5.10c)  and 
(3.10)  suggest  that  it  does  so  by  a  "small  amount"  that  vanishes  with  h  -  0. 

Since  the  initial  data  is  only  piecewise-  smooth  we  have  not  been  able  as  yet  to 
rigorize  these  arguments. 

We  remark  that  our  numerical  experiments  clearly  indicate  that  in  a  normal 
computational  situation  the  maximum  norm  of  the  numerical  solution  is  indeed 
uniformly  bounded.  We  feel  that  our  inability  to  prove  this  fact  stems  only  from 
lack  of  theoretical  tools  to  analyse  pointwise  regularity  of  the  numerical  solution. 


6.  The  nonlinear  case.  In  this  section  we  describe  an  approximate  solution 
v„(x,  f)  of  [5]  for  the  IVP  (2.2) 


(6.1)  v,  +  f(y)x  =  0,  v(x,0)  =  L(x;  v")  . 

This  approximate  solution  is  consistent  with  the  conservation  form  of  the  equation 

(6.1)  in  the  sense  that  the  cell-averaging  (2.3)  results  in  a  scheme  in  conservation 
form  i.e. 

(6.2)  v/+1  =  ~  f  v„(x,t)  dx  =  vj  -  \(fj+v2  ~  fj-in) 

where  the  numerical  flux  fj+m  is  consistent  with  /(«)  in  the  sense  of  (1.2c). 
Furthermore,  the  approximate  solution  operator  is  TVD 

(6.3)  7V(v„(-;  0)  s  7V(v„( ;  0))  =  7V(L(-;  v"))  for  0  as  r  S  t 

and  thus  by  the  reasoning  presented  in  section  2,  the  resulting  scheme  (6.2)  is 
non-osdllatory. 

We  turn  now  to  outline  the  derivation  of  this  approximate  solution.  To  sim¬ 
plify  our  presentation  we  ignore  entropy  considerations  and  refer  the  reader  to 
future  papers  for  details  of  appropriate  modifications.  We  assign  to  the  point 
xj+y 2  a  characteristic  speed  that  corresponds  to  the  Rankinc-Hugoniot  speed 
ciju.  1/2  of  the  two  neighbouring  cell-averages  vj  and  vj*+1 


(6.4a) 


aJ+V2 


'/(vjVi)  -  /(yfl 

V/+1  -  V/ 
a{yj) 


if  vf  *  v/+1 


and  denote  by  a(x)  the  piecewise-linear  interpolant  of  {fly +1/2},  i.e. 
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(6.4b)  a(x)  =  flji-i/2  +  («/+i/2  -  5)-i^)  •  (x  -  x,-i/2)/* 


for  X^_  y2  SlSi  X/4.V2  • 

The  approximate  solution  v„(x,t)  is  defined  by  specifying  its  constancy 
along  the  approximate  characteristics 

(6.5a)  x(0  =  x0  +  o(xq)  •  t 

i.e. 

(6.5b)  vH(x(t)t  t )  m  v„(x0,0)  =  L(xq;  vh)  . 

Using  (6.5a)  and  (6.4b)  to  express  jr0  in  terms  of  x  and  t 
6.5c)  xo(x,0  =  X)-m  +  h  [x  -  Xj-y&i) ]/[x>+1/2(f)  -  xj-viC 01 


for  xj-v&i)  <x<  xJ+v2(t) 


we  get  from  (6.5b)  that 
(6.6a)  v,,(x 


r,0  =  L  Xj_ 


-4-  h  .  - - - - .  y«] 

m  ^+l/2<0  ~  Xj-v2(t)  ’ 


for  x/.^r)  <  x  <  x^m(t) 


where 

(6.6b)  +1/2(0  =  *i+l/2  +  f  *  ®i+i/2  • 

Taking  cell-averages  of  the  approximate  solution  (6.6)  we  get  the  conservation 
form  (6.2) 

(6-7a)  v;+1  -  Vj  -  \(fj+  vi  ~  fj-yi) 


with  the  numerical  flux 


(6.7b) 


f(yf)  +  1/2  Sjf+i/2  (1  -  X  5/ -1/2)  ‘  Sj  if  ^j+1/2  s  0 
f(vJ+ 1)  -  1/2  1/2(1  +  X  «2;  +  3/2)  5y+i  if  aj+l/2  s  0  ’ 


where 

(6.7c)  5;  =  S]/[  1  +  X(a)+1/2  ~  «jf- 1/2)1  • 

Note  that  (6.7)  is  identical  to  (5.30)  in  the  constant  coefficient  case. 

We  turn  now  to  prove  that  the  scheme  (6.7)  is  uniformly  second-order  accu¬ 
rate  in  the  sense  of  (2.9).  Starting  with  the  exact  cell-averages  vj  -  uj  in  (6.7) 
this  amounts  to  showing  that 

(6*8a)  fj+1/2  =  fj+i/2(u)  +  ^(h2) 

with  a  sufficiently  smooth  coefficient  in  the  Oft2)  term;  here  fj4.  m  is  the 
numerical  flux  (6.7b)  computed  with  the  exact  cell-averages,  and  fj+w.(u)  is 
(2.1b).  We  shall  do  so  in  two  steps:  first  we  shall  show  that 

(6.8b)  fj+in(u)  =  y (/(£(*/+ 1/2>  «*"))  +  /(f-(xo(y/+  i/2>  t);  iT))J  x  Ofh2)  , 

where  *<](*./+ i/2>t)  is  (6.5c),  and  then  we  shall  verify  that 

(6.8c)  ^\f(L(x]+  2/2; «"))  +  f(L(xg[Xj+  m,  t);  «”))]  =  /y+1/2  +  0(*2)  . 

Special  attention  will  be  given  to  the  smoothness  of  the  Oft2)  coefficients. 

To  show  (6.8b)  we  start  by  using  the  trapezoidal  rule  to  approximate  the 
integral  in  (2.1b);  we  get 

(6.9a)  fj+iri(u)  =  y[fC«(x/+1/ 2,  f„))  +  /(u(xy+1/2,  tn  +  t))]  +  OQt2)  . 


The  smoothness  of  the  OQt2)  term  follows  from  that  of  /(it)  and  u(x,r).  Next 
we  observe  that  a(x)  in  (6.4b)  approximates  a(u(x,t„))  to  Oih2),  and  there¬ 
fore  we  can  use  the  approximate  characteristic  line  (6.5c)  to  trace 
«(*/+ V2>  *n  +  T)  to  u(xo(xj+V2>  T),  tn)  with  0(h2)  accuracy;  consequently 

(6.9b)  /(“(*/+ 1/2*  *n  +  t))  =  f(u(xo(Xj  + 1/2 »  t),  f„))  +  Oih2)  . 

Finally  we  obtain  (6.8b)  by  approximating  u(x,f„)  in  (6.9a)  and  (6.9b)  to 
Oih2)  by  L(x ;  «")  (see  (4.4)).  The  smoothness  of  the  ©(A2)  term  in  this 
approximation  is  due  to  (4.4c): 

SJ  -  A  •  «,(*.  »„)  +  £>(*3)  . 

(We  recall  that  the  degeneracy  to  first  order  accuracy  at  local  extrema  points  of 
some  "second-order  accurate"  TVD  schemes  is  due  to  lack  of  smoothness  there  of 
the  Oih2)  term  in  (2.7a)). 

We  turn  now  to  verify  (6.8c).  First  let  us  consider  the  case 

£(*/+i/2;  ^  + j  sj , 


£(*<)(*/+ 1/2> T);  «*")  -  uj  + 


Xa 


7+1/2 


1  +  X(fly+y2  “  <*j-V2) 


s? 


and  expand  the  LHS  of  (6.8c)  around  u].  We  get 

(6.10a)  fiuj)  +  j  a(57)(l  -  X  aj+yj)  Sj  +  -^'(^[(l  -  X  Sj 


+  (XJ^vj)2]^2  +  oih1)  =  //+w  +  y  (2X2a2  («JV«  +  0(A3) 


Similarly  in  the  case  aj+m  ^  0: 


£(*/+i/2;  «")  =  uj+i  -  j  5/+i,  Z.(jrQ(x>+La,  t);  5") 


_  .71  _ 

-  ty+1  ~ 


ka 


7+1/2 


2  1  +  \(aj+yi  -  ajf+i/j) 


■  f 


we  expand  the  LHS  of  (6.8c)  around  uf+i  to  get 
(6.10b)/(57+1)  -  —  a(uj+i)(l  +  kaj+y^Sj+i 


+  j  o'(»’*i)[(i  +  ^+m)2  +  (xs^w)2!  •  (£j+ 1)2  +  of*5) 

“  //-n/2  +  -j-  (2XJ“J  -  1)  •  a’  •  MV»  +  ^C*3)  ■ 

We  see  from  (6.10a)  and  (6.10b)  that  independently  of  the  sign  of  Ej+y 2,  the 
0(A2)  term  in  (6.8c)  is  the  same,  namely 

Y  (2X2«2  “  1)  *  '  (“*)2| /+1/2  • 

This  completes  the  proof  that  the  scheme  (6.7)  is  second-order  accurate  in  the 
sense  of  (2.9)  wherever  u(x,t)  is  smooth,  including  local  extrema  and  sonic 
if  =  0)  points. 

Remarks:  (1)  The  numerical  flux  (6.7b)  can  be  rewritten  as: 

(6.11)  f)+ 1/2  =  \  [/(v/)  +  /(v/+ 1)  ~  £/+l/2i  (V;>  1  -  vf) 


+  max(0,  fl/4.1/2)  •  (1  -  ^-1/2)  *  Sj 


-  min(0,  5} +1/2) 


(1  +  \aj+y 2)  • 
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(2)  Our  proof  that  the  scheme  (6.7)  is  non-oscillatory  is  based  on  the 
representation  of  (6.7)  as  the  cell-average  (6.2)  of  the  non-oscillatory  approxi¬ 
mate  solution  v„(x,f)  in  (6.6).,  To  ensure  that  vn(x,t)  remains  univalued  for 
Ostst  we  have  to  restrict  the  time-step  t  so  that  for  all  j 

(6.12a)  Xy+ 1/200  >  Xj-yii t)  . 

This  implies  the  condition 

(6.12b)  t  maXj(Sj-y2  -  5) +1/2)  =£  h  • 

On  the  other  hand,  to  derive  the  particular  form  of  the  numerical  flux  (6.7b)  we 
have  made  the  assumption 

(6.13a)  l*/+i/20O  ~  xj+uJi  <  h  for  all  j  , 

which  implies  the  condition. 

(6.13b)  t  maxyJSy+i^j  <  h  . 

The  two  time-step  restrictions  (6.12b)  and  (6.13b)  are  characteristic  to 
Godunov-type  schemes.  The  practice  of  reconciling  the  two  conditions  by 

t  maxj|fl)+1/2|  ^  ha 

is  completely  unnecessary:  The  scheme  (6.7),  as  the  original  Godunov  scheme, 
remains  non-oscillatory  (or  TVD  in  the  case  analysed  in  [10])  for  the  full  CFL- 
restriction  (6.13b).  The  reasoning  for  this  observation  is  as  follows:  The  approxi¬ 
mate  solution  (6.6)  under  the  CFL  restriction  (6.13b)  may  fail  to  be  univalued  in 
the  y'-th  cell  only  if  1/2  >  0  and  a}+i/2  <  0.  In  this  case  the  numerical  fluxes 
fj±\n  as  defined  by  cell-averaging  in  the  neighboring  cells  j  ±  1,  remain  (6.7b). 
Thus  the  only  thing  that  needs  to  be  changed  in  this  case  is  the  definition  of 
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vb(jc,t)  in  the  j-th  cell. 

(3)  We  observe  that  once  v„(jc.t)  is  defined  globally  in  (6.6)  there  is  no 
intrinsic  need  to  average  it  on  the  original  mesh.  We  may  average  it  on  different 
intervals  and  still  conclude  that  the  resulting  approximation  is  non-osdllatory  and 
conservative.  Furthermore,  the  construction  of  the  interpolant  Q,  the  approxi¬ 
mation  L  and  the  approximate  characteristic  field  a(x)  needed  to  define 
v„(x,  t),  does  not  depend  on  the  uniformity  of  the  mesh.  Therefore  the  scheme 
(6.7)  generalizes  immediately  to  non-uniform  moving  meshes.  Of  particular  com¬ 
putational  interest  are  the  self-adjusting  moving  grids  of  the  type  described  in 
[12],  which  make  it  possible  to  obtain  perfectly  resolved  shocks  and  contact 
discontinuities. 

(4)  We  note  that  since  the  approximate  solution  v„(x,t)  in  (6.6)  is  conser¬ 
vative,  it  is  possible  to  consider  an  associated  random-choice  method  obtained  by 
replacing  the  cell-averaging  in  (6.2)  by  a  sampling  with  respect  to  a  variable  that 
is  uniformly  distributed  in  the  cell,  i.e. 

*T‘  =  •’»(*/  + 

where  0"  is  uniformly  distributed  in  [-1/2, 1/2].  Following  the  reasoning  of  [7] 
it  is  clear  that  the  resulting  random-choice  method  is  non-osdUatory  and  that  its 
limits  are  weak  solutions  of  (1.1).  Although  the  random-choice  approach  has 
many  attractive  computational  features,  it  has  been  our  experience  that  in  many 
application  it  is  possible  to  accomplish  the  same  computational  goals  with  a  self- 
adjusting  moving  grid.  In  this  case  the  use  of  the  latter  is  preferable  as  it  offers 
gain  in  resolution  without  a  loss  in  pointwise  accuracy  that  is  assotiated  with  sam¬ 
pling. 
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7.  Numerical  Dhutration.  In  this  section  we  compare  the  new  uniformly 
second-order  non-osriUatory  scheme  of  this  paper  (to  be  refered  to  as  UN02)  to 
the  typical  second-order  TVD  scheme  (to  be  refered  to  as  TVD2).  Both  schemes 
can  be  written  in  the  form  (6.7),  i.c. 

(7.1a)  vj+1  -  vf  -  \(f>+1/2  ~  fj—1/2)  > 

/(vf)  +  1/2  ^(l  -  X5)_i^)SjV[l  +  X(a)+y2  “  5)— 1/2)1 

„  if  ^+1/2  0 

(71b)//.w  -  /(v?+l)  _  w  J/+1/2(X  +  xaJ+M>S;+1/tl  +  -  tyn/j)] 

if  flj+y2  ^  0 

(7.1c)  S;  =  m(5/^f); 

here  a^t-1/2  is  (6.4a)  and  m(x,y)  is  the  min  mod  function  (3.3).  Sf  are  dif¬ 
ferent  for  TVD2  and  UN02: 

(7.2)  TVD2:  Sf  ~dj±u2vn 

(7.3)  UN02:  5/  =  dJ±l/2vn  *  |  D,±1^v»  , 

where  and  D(+  ^  are  defined  in  (3.2). 

UN 02  and  TVD2  are  both  second-order  accurate  Godunov-type  schemes 
that  differ  only  in  the  reconstruction  step  (4.2a): 

(7.4a)  L(,xyu)  =  Uj  +  Sj(x  -  xj)Jh  for  [x  -  xj\  <  h/2  , 

where  the  slopes  of  the  lines  are  calculated  by  (7.3)  and  (7.2),  respectively. 
Therefore  we  start  our  comparison  on  the  approximation  level. 

In  Table  1  and  Fig.  1  we  present  approximations  to 
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u(x)  =  sinirx,  -lsjsl.  Wc  divide  [-1,1]  into  N  equal  intervals  and 
define 

(7.4b)  OSJSK. 

The  symbols  in  Fig.  1  denotes  values  of  uj  -  sin-n-xy  for  N  -  10  in  (7.4b).  In 
Fig.  la  we  show  the  piecewise-parabolic  interpolant  Q(x\u)  (see  section  3).  In 
Fig.  lb  we  show  the  piecewise-linear  approximation  Luno2(x;k)  which  is  (7.4a) 
with  (7.1c)  and  (7.3).  In  Fig.  lc  we  show  the  piecewise-linear  approximation 
lTVD^xju)  which  is  (7.4a)  with  (7.1c)  and  (7.2).  We  make  the  following  obser¬ 
vations  regarding  Fig.  1:  (i)  Q  is  a  better  approximation  than  IUN02;  LUN02  is  a 
better  approximation  than  L™2.  (ii)  7V(IUN02)  >  7V(u)  >  7V(LTVD2).  In 
table  1  we  quantify  the  first  observation;  we  list  the  L„-error  and  the  Lj -error  of 
these  approximations  to  simrx  for  a  refinement  sequence  of  N  -  10,20,40,80 
in  (7.4b).  Clearly  Q  is  an  0(h3)  approximation,  while  LUN02  and  LTVD2  arc 
Ofr2).  The  error  in  LUN02  is  about  a  LG  of  the  error  in  ltvdz. 

In  Table  2  and  Fig.  2  we  present  solutions  of  UN02  and  TVD2  for  the  con¬ 
stant  coefficient  case 

(7.6)  u,  +  ux  —  0,  m(x,0)  =  sin  irx,  -lsrsl 

with  periodic  boundary  conditions,  at  t  =  2  with  r/h  =  0.8.  Figs  3a  and  3b 
show  UN02  and  TVD2,  respectively,  with  N  -  20  in  (7.4b).  In  table  2  we  list 
the  I* -error  and  -error  for  a  refinement  sequence  with  N  =  20,40,80,160. 
Clearly  UN02  is  second-order  accurate  in  both  L»  and  L\,  while  TVD2  is 
second-order  accurate  in  L\  but  only  first  order  accurate  in  Lx.  Fig.  3b  demon¬ 
strates  that  the  degeneracy  to  first  order  accuracy  at  local  extrema  in  the  TVD 
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scheme  adversely  affects  the  accuracy  everywhere  (Because  the  scheme  is  TVD  it 
cannot  switch  abruptly  to  second-order  accuracy  as  this  would  introduce  osdlla- 
tious;  consequently  it  takes  quite  a  while  to  recover  the  second  order  accuracy). 

Next  we  approximate  the  discontinuous  function 

-xsin(3iTX2/2)  ,  -i  <  x  <  -1/3 
(7.7)  u(x)  =  |sin(2irx)|  ,  |x|  <  1/3 

2x  —  1  —  l/6sin(3'irx)  ,  1/3  <  x  <  1 

which  we  extend  to  have  period  2  outside  [- 1,1]. 

In  Fig  3  we  present  approximations  to  <j>(x),  using  N  -  20.  Fig  3a  shows 
Q(x,u),  Fig  3b  shows  L^^xjb),  and  Fig  3c  shows  LTVD2(x;u).  We  again 
observe  that  Q  is  better  approximation  then  LUN02,  while  L1^02  is  a  better 
approximation  then  LTVDZ. 

In  Fig  4  we  present  solutions  of  UNOZ  and  TVD2  for  the  constant  coeffi¬ 
cient  problem  (7.6),  initial  data  given  by  (7.2),  and  periodic  boundary  conditions. 
We  take  t  -  2  and  r/h  -  0.8.  Figs  4a  and  4b  show  UN02  and  TVD2  respec¬ 
tively  with  N  =  40.  Fig  4b  shows  the  damping  effect  that  the  TVD  scheme 
imposes  due  to  its  degeneracy  to  first  order  accuracy  at  local  extrema. 

In  Fig  5  we  solve  the  same  problems,  except  we  impose  boundary  conditions. 
At  x  =  -1  we  impose  the  given  function  (7.7)  evaluated  at  -1  -  t.  No  boun¬ 
dary  conditions  are  imposed  at  x  =  1.  We  implement  this  numerically  using 
UN02  and  TVD2  except  at  the  boundary  points.  There  we  are  in  general,  unable 
to  construct  non-osdllatory  piecewise  parabolic  interpolants  Q(x,u),  so  we  con¬ 
struct  the  only  possible  parabolic  interpolant  thru  x„xl+i  and  the  point  to  either 
the  left  or  right  which  lies  in  the  region.  The  analogous  procedure  is  carried  out 


at  the  reconstruction  stage.  Figs  5a  and  5b  again  show  the  results  at  t  =  2  with 
r/h  =  0.8. 


Hie  possible  introduction  of  oscillations  through  the  boundary  conditions 
does  not  seem  to  have  degraded  the  performance  of  either  scheme  (in  fact  the 
opposite  is  observed).  Again  the  TVD2  scheme  shows  a  damping  effect. 

hi  Table  3  and  Fig.  6  we  present  results  for  Burgers’  equation 

(7.8)  u,  +  uitx  -  0,  u(x,0)  =  a  +  sin  ir(x  +  0),  -lsxs  1 

with  periodic  boundary  conditions  and  r/h  (1  +  |a|)  =  0.5.  The  solution  to  (7.7) 
is  smooth  for  t  <  1/ir;  at  t  —  Vrt  it  develops  shocks.  In  Table  3a  we  list  the 
Loo-error  and  Lj-error  of  UN02  and  TVD2  at  /  =  0.15  for  a  =  0  =  0  in 
(7.7).  This  table  shows  the  same  asymptotic  behaviour  as  Table  2,  except  that 
because  of  the  large  gradients  it  shows  for  a  smaller  h. 

In  Figs  6a  and  6b  we  show  results  of  UN02  and  TVD2  for  (7.8)  with  a  =  2 
and  0  =  1  at  t  =  0.35  with  N  —  20.  In  this  case  the  solution  to  (7.8) 
develops  a  shock  moving  with  speed  2  beginning  at  time  t  -  1/tt  ~  0.318. 

In  Table  3b  and  Figs  6c  and  6d  we  repeat  the  previous  calculations  for  the 
schemes  (2.13): 

(7  *>)  ^  ~  ^</j+V2  ~  fi-m) 


(7'9b)  fj*ui  -  /(*„(*/+ 1/J,  x/2))  =  /(i(x0(ly+ V2,  T/2),v")> 

where  xd(xj+\f2,  r/2)  is  (6.5c),  i.e. 


(7 .9c)fj+i/2  - 


/ 

/ 


L  +  I  .  1  ~  +  &i-inyi  n 

l  2  1  +  \(&j+v2  -  2)/2  J 

(„  _  1  +  H&J+V2  +  . 

1  2  i  +  X(<Sy+ ~  bj+viyi 


if  &J+V2  a  0 
^  &j+  i/2  s  0 


As  we  have  remarked  in  section  2,  v"+1  in  (7.9a)  is  not  a  cell-average  of 
v„(x,t),  but  only  an  approximation  to  it.  Therefore  it  is  not  necessary  to  take 
&j+u2  in  (7.8c)  to  be  (6.4a).  We  choose  so  that  (7.9c)  is  continuous  at 

&j+V7  -  0: 

(7.9d)dy+v2  = 

We  denote  the  schemes  (7.9)  with  S]  defined  by  (7.1c)  and  either  (7.2)  or 
(7.3)  by  FVD2  and  FN02,  respectively.  We  note  that  (7.9)  is  identical  to  (7.1)  in 
the  constant  coefficient  case,  and  consequently  FVD2  and  FN02  are  nonoscilla- 
tory  in  the  constant  coefficient  case.  Figs  $c  and  6d  show  that  FN02  and  FVD2 
are  also  non  oscillatory  in  the  case  (7.8).  Furthermore,  Table  3b  shows  that 
FN02  is  much  more  accurate  than  UN02  (FVD2  is  about  the  same  as  TVD2). 

hi  all  previous  examples  we  have  presented  pointwise  calculations;  namely, 
we  have  initialized  the  numerical  solution  by  taking  vj  to  be  the  value  of  the  ini¬ 
tial  data  at  xj,  and  we  have  considered  v"  to  be  an  approximation  to  u{xjytn). 
(Surely  this  is  an  acceptable  practice  for  second  order  accurate  schemes.)  In  Table 
3c  we  repeat  the  calculation  for  UN02  in  Table  3a,  but  now  in  a  sense  of  cell- 
averages  and  denote  it  by  AN02.  Now  we  initialize  UN02  for  (7,8)  with 
a  =  3  -  0  by  cell-averages  of  the  initial  data,  i.e. 


/I 


<»Ai  -  -  M  +  I SP 


(»Ai  -  -  0*  +  \sj) 


(7.10a)vj>  =  [cos(irxy ..  1/2)  -  cot^itXj— 1/2)1  =  sin(uxy)  , 

and  regard  v"  to  represent  cell-averages  of  u  (*,/„).  To  obtain  a  pointwise 
approximation  to  u(x,tn)  we  first  compute  point  values  VJ+V2  of  its  indefinite 

integral  =  f  /*”u(y,fII)  <fy  by 

x0 
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(7.10b) 


«-<0 


Next  we  obtain  a  global  piecewise- linear  approximation  v(x,tn)  to  u(x,tn )  by 
(7.10c)  v(x,r„)  =  ~  C(x;V") 

where  Q  is  the  piecewise-parabolic  interpolant  of  section  3.  Finally  we  get 
(7.i0d)  v(xj,tn)  -  £  Q(xfy)  =  i  (v;+w  -  v;-w)  =  v/ . 

Thus  the  only  difference  between  AN02  in  Table  3c,  and  UN02  in  Table  3a 
is  the  initialization  (7.10a),  which  itself  differs  only  slightly  from  the  mesh  values 
of  the  initial  data  (since  sin(irA/2)/(ir/i/2)  =  1  -  l/6(vh/2)2  +  0(A4)). 

We  remark  that  cell-averages  do  play  a  significant  role  when  the  initial  data 
is  discontinuous  (since  they  provide  information  about  the  location  of  the  discon¬ 
tinuity)  and  in  higher-order  Godunov-type  schemes;  this  will  be  described  else¬ 
where. 

We  turn  now  to  present  calculations  with  a  formal  extension  of  UN02  and 
TVD2  for  systems  of  conservation  laws.  We  consider  a  Riemann  problem  for  the 
Euler  equations  of  gas  dynamics 

(7.11a)  u,  +  f{u)x  =  0,  u(x,0)  = 

(7.11b)  u  =  (p,mJE)T,  f(u)  =  (m,m2/p  -I-  P,  m(E  +  P)/ p)T  , 

(7.11c)  P  =  h  -  1  )(E  -  i  m2/p)  . 

Here  p,  m,  E  and  P  are  the  density,  momentum,  total  energy  and  pressure, 


x  <  0 

w 

x  >  0 
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except  with  given  data  Figure  5b.  Same  as  4b  except  with  given  data 

and  outflow  at  x  =  1 .  at  x  =  -1  and  outflow  at  x  =  1 


1.26 


given  (7.7)  with 


Figure  2.  Numerical  solutions  of  u+  +  u  =  0,  u(x,0)  -  simrx  at  t  =  2,  with 


PIECEWISE  LINERR  TVD  PIECEWISE  LINERR 


U_  ra 


PIECEWISE  PARABOLIC 


Figure  la 

igure  1 .  Approximations  of  u  =  simrx,  -1  <  x  <  1,  with  N 
)  Piecewise-Parabolic  interpolant  Q(X,u). 


Table  3a.  Numerical  solutions  of  u,  +  uux  —  0,  u(x,0)  =  sinirr,  at  t  -  0.15  and  t/h 
UN02  and  TVD2 


L  » -ERROR 

Lr ERROR 

UN02 

TVD2 

UN02 

TVD2 

40 1 5  712  x  10-3  1054  X  10  l 3  034x  10  5  051  x  10 


1.552  x  10"3  4.422  x  10"3  |  7.771  x  10~4  1.340  x  10“3 


1 160  |  3.985  x  10"4 
Table  3b.  Same  as  Table  3a  for  FN02  and  FVD2 


Loo-ERROR 

Lj -ERROR 

FN02 

FVD2 

FN02 

FVD2 

1.959  x  10" 3  1.054  x  10"2 


4.424  x  10"3 


1.837  x  10“3 


3.835  x  10-3 


1.072  X  10"3 


Table  3c.  Same  as  table  3a  for  AN02. 


L  ao-ERROR  |  Lj -ERROR 


2.249  x  10"2  I  1.221  x  10" 2 


6.623  x  10"3  I  3.243  x  10"3 


1.781  x  10"3  I  8.259  x  10"4 


160  I  4.597  x  10"4  I  2.079  x  10 


Table  1:  Approximations  to  u(x)  -  sin(irac),  -lsxsl,  with  periodic  boundary  conditions. 


Likewise  5*  denotes  the  component  of  the  vector  of  slopes  in  the  4-th  charac¬ 
teristic  field,  and  is  defined  as  follows: 

(7.13e)  Sf  =  m(S»  +  HSj+m  -  ^-w)] ; 

m(x,y )  is  the  min  mod  function  (3.3).  S±  j  are  different  for  VD2  and  N02: 

(7.14)  TVD2: 5^  =  d}±y2w 

(7.15)  UN02:  Sk±J  =  df+mw  =p  y 

£>*±1/2*  =  n»(df+3/2w  -  df+y2w,  df+y2w  -  . 

In  Figs  7.8  and  7.9  we  show  numerical  solutions  of  UN02  and  TVD2,  respec¬ 
tively,  for  the  Riemann  problem  (7.1b)  with 

UL  =  (l,0,2.5)r,  UR  =  (0.125,0,0.25)  . 

These  figures  demonstrate  that  the  formal  extension  to  systems  is  nonosdllatory 
in  this  case.  Since  the  solution  to  the  Riemann  problem  is  just  constant  states 
seperated  by  waves  we  do  not  get  to  see  here  the  extra  resolution  power  of  UN02, 
except  that  its  numerical  solution  is  somewhat  '’crisper'’  than  that  of  TVD2.  In 
this  calculation  we  have  not  employed  any  artificial  compression  in  the  linearly 
degenerate  field  and  therefore  the  contact  discontinuity  smears  like  n173,  as 
expected.  The  interested  reader  is  referred  to  [4],  [5]  and  [10]  for  a  detailed 
description  of  such  compression  techniques,  as  well  as  for  details  of  entropy 
enforcement  mechanisms. 


respectively;  we  take  y  =  1.4. 

In  the  following  we  apply  the  extension  technique  of  [3]  to  UN02  and  TVD2. 
The  idea  is  to  extend  UN02  and  TVD2  to  systems  in  such  a  way  that  will  be 
identical  to  (7.1)  in  the  scalar  case,  and  will  decouple  into  (5.3)  for  each  of  the 
characteristic  variables  in  the  constant  coefficient  system  case.  To  accomplish 
that  we  use  Roe’s  averaging  for  (7.11)  (see  [13]) 

(7.12a)  Vj+y2  -  V(yf,  v/+1) 

for  which 

(7.12b)  fiyf+i)  -  f(yj)  =  A(v/+y2)(v;+1  -  vjO,  A(u)  =  df/du  , 

and  define  local  characteristic  variables  with  respect  to  the  right-eigenvector  sys¬ 
tem  of  A(yj+1/2).  We  extend  (6.11)  to  systems  as  follows: 

(7.13a)  v;+1  -  vj  -  Xtfj+M  -  fj-vd 

(7. 13b)  W  =  \  U-f)  +  flyf* ,)  -  cf+  m  1 m 

(7.13c )'cf+V2.  "  \3}+vi]4+i nw  “  max(0,fl/+y2)(l  - 

+  min(0,5^+ !/2)(l  +  \aj+y2)Sj+i  • 

Here  aj+y 2  is  the  t-th  eigenvalue  of  A(vj+y 2)  corresponding  to  and 

df+  yzw  denotes  the  component  of  dj+  yjy  =  v"+1  -  v"  in  the  *-th  characteris¬ 
tic  field,  i.e. 

3  . 

dj+V2v  ~  2  (dJ+V 2w)R}+V2  • 


(7.13d) 
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